*** to be run after 1_Table_2\7_estimation_sample_and_variables.do
cap log close	
	
log  using  ".../AA_8_ben_08Sep2020.smcl",replace

version 13	
use .../sample_estimation_1_03Sept2020.dta, clear
		

***** generating binary variable if individual has received any unemployment benefit during the year
gen unemp= unemp_ben > 0 & unemp_ben != .

gen se_unemp  =unemp*sedum*(1-lcens_spell)
gen self_unemp=unemp*sedum*(lcens_spell)
gen e_unemp   =unemp*(1-sedum)*(1-lcens_spell)
gen elf_unemp =unemp*(1-sedum)*(lcens_spell)

*** 1 ***	
	cap program drop prog2_lnf
	set type double					

	program define prog2_lnf
*       version 8.2

	args todo b lnf

	tempvar ISE IU ISELF IULF logitp2 logitp3 logitp4 b1 b2 b3 b4   a1 a2 a3 a4 logitp5 logitp6 logitp7 logitp8  lnfie last lE1 lE2 lE3   lE1LF lE2LF lE2LF ///
	        lU1  lU2 lU3 lU1LF lU2LF lU3LF l1 ll1 l2 ll2 l3 ll3 l4 ll4 l5 ll5 l6 ll6 l7 ll7 l8 ll8  mlSE1 mlSE2  mlSE3  mlE1 mlE2 mlE3 mlSELF1  mlSELF2 mlSELF3  mlELF1 mlELF2 mlELF3 ///
	        h1SE h2SE h3SE  sum1SE sum2SE  sum3SE  ST1SE ST2SE ST3SE lE3LF h1SELF h2SELF  h3SELF  sum1SELF sum2SELF  sum3SELF  ST1SELF ST2SELF ST3SELF ///
		h1U h2U  h3U  sum1U sum2U  sum3U  ST1U ST2U ST3U h1ULF h2ULF   h3ULF  sum1ULF sum2ULF  sum3ULF  ST1ULF ST2ULF ST3ULF lastSE lastU lastSELF lastULF lnfi
		
	tempname p2 p3 p4 p1
	
        *defining equations
    mleval `ISE'     = `b', eq(1)
	mleval `IU'      = `b', eq(2)
    mleval `ISELF'   = `b', eq(3)
	mleval `IULF'    = `b', eq(4)
	mleval `a1'       = `b', scalar eq(5)
	mleval `a2'       = `b', scalar eq(6)
	mleval `b1'       = `b', scalar eq(7)
	mleval `b2'       = `b', scalar eq(8)
	mleval `logitp2' = `b', scalar  eq(9) //5



	

	quietly {
			* to make sure probabilities constrained between 0 and 1
			* for identification, p11==1-p21-p22 and m11==0
			scalar `p2' = exp(`logitp2')/(1+exp(`logitp2'))
			scalar `p1' = 1             /(1+exp(`logitp2'))

			***** for SELFEMPLOYED spells
			
			* hazard functions for the three types
			by $S_E_id $spell: gen double `h1SE' = 1-exp(-exp(`ISE' + `a1'))	if $SE_E ==1 & $LF == 0
			by $S_E_id $spell: gen double `h2SE' = 1-exp(-exp(`ISE' ))			if $SE_E ==1 & $LF == 0
			
                        * first part of survival function 
			by $S_E_id $spell: gen double `sum1SE' = sum(exp(`ISE' + `a1'))		if $SE_E ==1 & $LF == 0
			by $S_E_id $spell: gen double `sum2SE' = sum(exp(`ISE' ))			if $SE_E ==1 & $LF == 0
		
			
			* cumulative product of survival functions
			by $S_E_id $spell: gen double `ST1SE' = exp(-`sum1SE'[_N])			if _n==_N & $SE_E ==1 & $LF == 0
			by $S_E_id $spell: gen double `ST2SE' = exp(-`sum2SE'[_N])			if _n==_N & $SE_E ==1 & $LF == 0
			

			* identify last obs of the spell
			by $S_E_id $spell: gen byte `lastSE' = (_n==_N) if $SE_E ==1 & $LF == 0
			replace `lastSE'=0 if `lastSE'==.

			* likelihood functions of the spells for each type
			by $S_E_id $spell: gen double `lE1' =  `ST1SE' * ((`h1SE'/(1-`h1SE'))^$event ) if _n==_N & $SE_E ==1 & $LF == 0
			by $S_E_id $spell: gen double `lE2' =  `ST2SE' * ((`h2SE'/(1-`h2SE'))^$event ) if _n==_N & $SE_E ==1 & $LF == 0
			
			
			replace `lE1'=0 if `lE1'==.
			replace `lE2'=0 if `lE2'==.
			

			******** for EMPLOYED 
			by $S_E_id $spell: gen double `h1U' = 1-exp(-exp(`IU' + `b1' )) if $SE_E ==0 & $LF == 0
			by $S_E_id $spell: gen double `h2U' = 1-exp(-exp(`IU' )) if $SE_E ==0 & $LF == 0
			
			
			by $S_E_id $spell: gen double `sum1U' = sum(exp(`IU' + `b1' )) if $SE_E ==0 & $LF == 0
			by $S_E_id $spell: gen double `sum2U' = sum(exp(`IU' )) if $SE_E ==0 & $LF == 0
			

			by $S_E_id $spell: gen double `ST1U' = exp(-`sum1U'[_N]) if _n==_N & $SE_E ==0 & $LF == 0
			by $S_E_id $spell: gen double `ST2U' = exp(-`sum2U'[_N]) if _n==_N & $SE_E ==0 & $LF == 0
			
			
			by $S_E_id $spell: gen byte `lastU' = (_n==_N) if $SE_E ==0 & $LF == 0
			replace `lastU'=0 if `lastU'==.

			by $S_E_id $spell: gen double `lU1' =  `ST1U' * ((`h1U'/(1-`h1U'))^$event ) if _n==_N & $SE_E ==0 & $LF == 0
			by $S_E_id $spell: gen double `lU2' =  `ST2U' * ((`h2U'/(1-`h2U'))^$event ) if _n==_N & $SE_E ==0 & $LF == 0
                        
			               
			replace `lU1'=0 if `lU1'==.
			replace `lU2'=0 if `lU2'==.
			
			
			***** for SELFEMPLOYED left-censored spells
			
			* hazard functions for the three types
			by $S_E_id $spell: gen double `h1SELF' = 1-exp(-exp(`ISELF' + `a2'))   if $SE_E ==1 & $LF == 1
			by $S_E_id $spell: gen double `h2SELF' = 1-exp(-exp(`ISELF'))   if $SE_E ==1 & $LF == 1
			
			
                        * first part of survival function
			by $S_E_id $spell: gen double `sum1SELF' = sum(exp(`ISELF' + `a2'))  if $SE_E ==1 & $LF == 1
			by $S_E_id $spell: gen double `sum2SELF' = sum(exp(`ISELF' ))    if $SE_E ==1 & $LF == 1
			
			
			* cumulative product of survival functions
			by $S_E_id $spell: gen double `ST1SELF' = exp(-`sum1SELF'[_N]) if _n==_N & $SE_E ==1 & $LF == 1
			by $S_E_id $spell: gen double `ST2SELF' = exp(-`sum2SELF'[_N]) if _n==_N & $SE_E ==1 & $LF == 1
			
			
			* identify last obs of the spell
			by $S_E_id $spell: gen byte `lastSELF' = (_n==_N) if $SE_E ==1 & $LF == 1
			replace `lastSELF'=0 if `lastSELF'==.

			* likelihood functions of the spells for each type
			by $S_E_id $spell: gen double `lE1LF' =  `ST1SELF' * ((`h1SELF'/(1-`h1SELF'))^$event ) if _n==_N & $SE_E ==1 & $LF == 1
			by $S_E_id $spell: gen double `lE2LF' =  `ST2SELF' * ((`h2SELF'/(1-`h2SELF'))^$event ) if _n==_N & $SE_E ==1 & $LF == 1
			
           
			replace `lE1LF'=0 if `lE1LF'==.
			replace `lE2LF'=0 if `lE2LF'==.
			
			

			******** for EMPLOYED left-censored spells
			by $S_E_id $spell: gen double `h1ULF' = 1-exp(-exp(`IULF' + `b2')) if $SE_E ==0 & $LF == 1
			by $S_E_id $spell: gen double `h2ULF' = 1-exp(-exp(`IULF' )) if $SE_E ==0 & $LF == 1
			

			by $S_E_id $spell: gen double `sum1ULF' = sum(exp(`IULF' + `b2')) if $SE_E ==0 & $LF == 1
			by $S_E_id $spell: gen double `sum2ULF' = sum(exp(`IULF' )) if $SE_E ==0 & $LF == 1
			
			
			by $S_E_id $spell: gen double `ST1ULF' = exp(-`sum1ULF'[_N]) if _n==_N & $SE_E ==0 & $LF == 1
			by $S_E_id $spell: gen double `ST2ULF' = exp(-`sum2ULF'[_N]) if _n==_N & $SE_E ==0 & $LF == 1
			
			
			by $S_E_id $spell: gen byte `lastULF' = (_n==_N) if $SE_E ==0 & $LF == 1
			replace `lastULF'=0 if `lastULF'==.

			by $S_E_id $spell: gen double `lU1LF' =  `ST1ULF' * ((`h1ULF'/(1-`h1ULF'))^$event ) if _n==_N & $SE_E ==0 & $LF == 1
			by $S_E_id $spell: gen double `lU2LF' =  `ST2ULF' * ((`h2ULF'/(1-`h2ULF'))^$event ) if _n==_N & $SE_E ==0 & $LF == 1
			
                                       
			replace `lU1LF'=0 if `lU1LF'==.
			replace `lU2LF'=0 if `lU2LF'==.
			
	    
				// individual log likelyhood for each type
			by $S_E_id: gen double `l1'=`lE1' +`lU1'+`lU1LF'+`lE1LF'
			by $S_E_id: gen double `ll1'=`p1'*exp(sum(log(`l1')))
			by $S_E_id: gen double `l2'=`lU2'+`lE2' +`lU2LF'+`lE2LF'
			by $S_E_id: gen double `ll2'=`p2'*exp(sum(log(`l2')))




			by $S_E_id: gen byte `last' = (_n==_N) 
	
			* individual contribution to the log likelihood 
			gen double `lnfie' = cond( !`last'  , 0, $weight * ln(( `ll1') + ( `ll2')   ))
				
					
			mlsum `lnf' = `lnfie'

	}
end

	
	
	

set type double		


tempvar mysamp 
		tempname b b01 b02 b1 b2 b03 b04 b3 b4 lnf V  a1 a2 a3  a4 b4 masshE masslE masshSE masslSE ///
		       logitp2 logitp3 logitp4  masslSELF masslELF  masshSELF masshELF

global varse " se_netincdiff se_convexity se_female_1 se_ageb se_edu_2-se_edu_3  se_yr_2-se_yr_4 se_reg_2-se_reg_5 se_unemp"
global vare  " e_netincdiff e_convexity e_female_1 e_ageb e_edu_2-e_edu_3  e_yr_2-e_yr_4 e_reg_2-e_reg_5 e_unemp"
global varself " self_netincdiff self_convexity  self_female_1 self_ageb self_edu_2-self_edu_3 self_yr_2-self_yr_4 self_reg_2-self_reg_5 self_unemp"
global varelf  "elf_netincdiff elf_convexity elf_female_1 elf_ageb elf_edu_2-elf_edu_3  elf_yr_2-elf_yr_4 elf_reg_2-elf_reg_5 elf_unemp"
global dumse "  lnt_se"
global dume "    lnt_e"
global dumself " lnt_self "
global dumelf " lnt_elf"

 
set more off	
// initial values for ll function taken from clogclog		
glm durvar  $varse  $dumse if sedum==1 & lcens_spell ==0 ,  f(b) l(cloglog) 
matrix `b01' = e(b)
matrix coleq `b01' = hazardE
	
glm durvar   $vare  $dume  if sedum==0 & lcens_spell ==0   ,  f(b) l(cloglog) 	
matrix `b02' = e(b)
matrix coleq `b02' = hazardU

glm durvar  $varself  $dumself  if sedum==1 & lcens_spell ==1 ,  f(b) l(cloglog) 
matrix `b03' = e(b)

matrix coleq `b03' = hazardELF
	
glm durvar   $varelf  $dumelf  if sedum==0 & lcens_spell ==1 ,  f(b) l(cloglog) 	
matrix `b04' = e(b)

matrix coleq `b04' = hazardULF

local LL1 = e(ll)
// create matrix with initial betas from cloglog and initial mass point location and probability (given)

matrix `a1' = (-0.1 )
matrix colnames `a1' = a1:_cons

matrix `a2' = (-1.3)
matrix colnames `a2' = a2:_cons

matrix `a3' = (-3.2)
matrix colnames `a3' = a3:_cons

matrix `b1' = (-0.2 )
matrix colnames `b1' = b1:_cons

matrix `b2' = (-1.2 )
matrix colnames `b2' = b2:_cons
matrix `b3' = (-1.5 )
matrix colnames `b3' = b3:_cons


matrix `logitp2' = (0.5)
matrix colnames `logitp2' = logitp2:_cons
matrix `logitp3' = (-2.1)
matrix colnames `logitp3' = logitp3:_cons
matrix `logitp4' = (-1.1)
matrix colnames `logitp4' = logitp4:_cons


matrix `b1' = `b01',`b02',`b03',`b04',`a1',`a2',`b1',`b2',`logitp2'


matrix list `b1'

// define macro that identify id variable
global S_E_id "id"
// define macro that identify se or E spell
global SE_E "sedum"
// define macro that identy the spell
global spell "a_spell"
// define macro that identify event
global event "durvar" 
 // define macro that identify left_censored spells
global LF "lcens_spell" 
 
// define weight
global weight "pweight"
 
 
sort id a_spell year
 
*cap log close
*log  using /ssb/ovibos/a2/ekstern/apa/prog/apa/survival/moredefition/durmod_`gruppo'.log,replace

//ml model, d0, no gradient or Hessian provided
ml model d0 prog2_lnf (hazardE: durvar = $varse  $dumse ) ///
                       (hazardU: durvar =  $vare  $dume  ) ///
		       (hazardELF: durvar = $varself $dumself   ) ///
		       (hazardULF: durvar = $varelf $dumelf  ) ///
					    (a1: ) (a2: ) (b1: ) (b2: )    ///
					    (logitp2: )   ///
                                , waldtest(0) /*init(`b1')  */ 


ml init hazardE:se_netincdiff    		=  -.42 
ml init hazardE:se_convexity    		=   .04 
ml init hazardE:se_female_1    			=  -.02 
ml init hazardE:se_ageb    				=  -.01 
ml init hazardE:se_edu_2    			=   .00 
ml init hazardE:se_edu_3    			=   .21 
ml init hazardE:se_yr_2    				=  -.09 
ml init hazardE:se_yr_3    				=   .02 
ml init hazardE:se_yr_4    				=  -.22 
ml init hazardE:se_reg_2    			=   .03 
ml init hazardE:se_reg_3    			=   .01 
ml init hazardE:se_reg_4    			=   .03 
ml init hazardE:se_reg_5    			=   .09 
ml init hazardE:se_unemp    			=   .26
ml init hazardE:lnt_se    				=  -.51 
ml init hazardE:_cons    				= -1.13 
ml init hazardU:e_netincdiff    		=   1.68
ml init hazardU:e_convexity    			=   -.24
ml init hazardU:e_female_1    			=    .60
ml init hazardU:e_ageb    				=    .02
ml init hazardU:e_edu_2    				=    .11
ml init hazardU:e_edu_3    				=    .13
ml init hazardU:e_yr_2    				=    .24
ml init hazardU:e_yr_3    				=    .30
ml init hazardU:e_yr_4    				=   -.14
ml init hazardU:e_reg_2    				=   -.22
ml init hazardU:e_reg_3    				=   -.32
ml init hazardU:e_reg_4    				=   -.41
ml init hazardU:e_reg_5    				=   -.47
ml init hazardU:e_unemp    				=   .44	
ml init hazardU:lnt_e    				=   -.49
ml init hazardU:_cons    				=  -3.10
ml init hazardELF:self_netincdiff    	=  -.72
ml init hazardELF:self_convexity    	=  -.01
ml init hazardELF:self_female_1    		=   .19
ml init hazardELF:self_ageb    			=  -.03
ml init hazardELF:self_edu_2    		=  -.13
ml init hazardELF:self_edu_3    		=   .05
ml init hazardELF:self_yr_2    			=  -.16
ml init hazardELF:self_yr_3    			=  -.31
ml init hazardELF:self_yr_4    			=  -.73
ml init hazardELF:self_reg_2    		=   .03
ml init hazardELF:self_reg_3    		=  -.04
ml init hazardELF:self_reg_4    		=   .18
ml init hazardELF:self_reg_5    		=   .13
ml init hazardELF:self_unemp    		=   .82
ml init hazardELF:lnt_self    			=  -.01
ml init hazardELF:_cons    				=  -.89
ml init hazardULF:elf_netincdiff    	=   1.75
ml init hazardULF:elf_convexity    		=   -.16
ml init hazardULF:elf_female_1    		=    .77
ml init hazardULF:elf_ageb    			=   -.04
ml init hazardULF:elf_edu_2    			=   -.00
ml init hazardULF:elf_edu_3    			=    .10
ml init hazardULF:elf_yr_2    			=    .15
ml init hazardULF:elf_yr_3    			=    .03
ml init hazardULF:elf_yr_4    			=   -.74
ml init hazardULF:elf_reg_2    			=   -.29
ml init hazardULF:elf_reg_3    			=   -.48
ml init hazardULF:elf_reg_4    			=   -.55
ml init hazardULF:elf_reg_5    			=   -.53
ml init hazardULF:elf_unemp    			=    .92
ml init hazardULF:lnt_elf    			=   -.23
ml init hazardULF:_cons    				=  -1.92
ml init a1:_cons    					=  -0.53
ml init a2:_cons    					=  -0.33 
ml init b1:_cons    					=  -0.04   
ml init b2:_cons    					=  -0.83
ml init logitp2:_cons    				=  -0.41
					
ml maximize, trace  search(off) difficult gradient

eststo m


eststo r1lr: nlcom (p1: 1/(1+exp([logitp2]_cons))) ///
		   (p2: exp([logitp2]_cons)/(1+exp([logitp2]_cons))), post


log close
